## To run SRE model over a range of range-edge gradients, dispersal costs, and Repro rates

source(file="/home1/30/jc227089/evo-dispersal/SRE/Stable range functions.R", echo=F)
setwd("/home1/30/jc227089/SRE/gvar_adap/ran_outs")
load("/home/30/jc227089/SRE/equil/ran_outs/SRE-equil_summary.RData")
### Functions ###

# to collect data
collect<-function(popmatrix, spX, spY, gen, hab.mat, ...){ #edge x is last sustainable population column
	if (length(popmatrix[,1])==0) return(NULL)	
	popsize<-tapply(popmatrix[,"P"], popmatrix[,"X"], length) #popsize over columns
	mean.disp<-tapply(popmatrix[,"P"], popmatrix[,"X"], mean) #mean trait value over columns
	mean.h<-tapply(popmatrix[,"H"], popmatrix[,"X"], mean)
	temp<-hab.mat[popmatrix[,"X"], 1]
	malad<-abs(popmatrix[,"H"]-temp)
	malad<-tapply(malad, popmatrix[,"X"], mean) #mean maladaptation over columns
	gen<-rep(gen, length(popsize))
	out<-list(c(...), cbind(gen, popsize, mean.disp, mean.h, malad))
	out
}

gvar_adap<-function(n=20*40*20, ngens.init=250, ngens=1250, maxX=40, spY=20, end.samp=20, steps, K=40, lambda, d.cost, mutn=0.05, plot=F, dval=99){
	pop<-init.inds(n, maxX, spY, hvar=T, dval)
	if (dval<=1) dmut<-FALSE
	else dmut<-TRUE
	hab<-init.hab(maxX, spY, 0)
	nbrs<-neighbours.init(maxX, spY) #reveal all neighbours
	data<-list(par=c(), space=c()) # to take last 50 gens
	for (i in 1:ngens.init){
		pop<-repro.disp(pop, hab, maxX, spY, K, lambda, d.cost, nbrs, mutn, adap=T, dmut=dmut)
		if (plot==T) plotter(pop, maxX, spY, K)
	}
	pop<-subset(pop, pop[,"X"]<=maxX/2) #remove half the population
	hab<-init.hab(maxX, spY, steps) #drop habitat gradient over the population
	for (i in (ngens.init+1):ngens){
		pop<-repro.disp(pop, hab, maxX, spY, K, lambda, d.cost, nbrs, mutn, adap=T, dmut=dmut)
		if (plot==T) plotter(pop, maxX, spY, K)
		if (ngens-i<=end.samp) {
			temp<-collect(pop, maxX, spY, i, hab, steps, lambda, d.cost, dval)
			data[[1]]<-rbind(data[[1]], temp[[1]])
			data[[2]]<-rbind(data[[2]], temp[[2]])
		}
	}
	colnames(data[[1]])<-c("steps", "lambda", "d.cost", "dval")
	xloc<-as.numeric(row.names(data[[2]]))
	data[[2]]<-cbind(xloc, data[[2]])
	data
}

### initialise ###
args=(commandArgs(TRUE))

#evaluate the arguments
# input arguments will be number of replicates runs (rr), file.ID, lambda (lv), steps (stps) and dispersal cost (dc)
for(i in 1:length(args)) {
	 eval(parse(text=args[[i]]))
}

pars<-c(lv, stps, dc, dvl)
sim.result<-c()
for (rep in 1:rr){
	out<-gvar_adap(steps=stps, lambda=lv, d.cost=dc, dval=dvl)
	out<-c(list(pars), out)
	sim.result<-cbind(sim.result, out)
}
save(sim.result, file=paste("gvarout", file.ID, ".RData", sep=""))